Coarse-graining the dynamics of coupled oscillators 
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We present an equation-free computational approach to the study of the coarse-grained dynamics 
of finite assemblies of non-identical coupled oscillators at and near full synchronization. We use 
coarse-grained observables which account for the (rapidly developing) correlations between phase 
angles and oscillator natural frequencies. Exploiting short bursts of appropriately initialized detailed 
simulations, we circumvent the derivation of closures for the long-term dynamics of the assembly 
statistics. 
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Since Winfree's pioneering work in 1960's 0], cou- 
pled oscillator models have been investigated exten- 
sively. Some exact results on the collective dynamics for 
an infinite number of coupled oscillators (the so-called 
continuum-limit) EJ have shed light on synchro- 
nization phenomena in biological IH l a IS EI B , chemi- 
cal HUHl! and physical systems [ll|,[l2|]- However, even 
in this ideal limit, some basic questions including global, 
quantitative stability of asymptotic states, still remain 
open 0, 0, . Many real- world systems consist 
of a large, finite number of non-identical entities, where 
statistical techniques for the continuum-limit are not di- 
rectly applicable. Exploring and understanding the dy- 
namics of such finite oscillator assemblies is an important 
topic (e.g., see Ref. |3). 

We present a computer-assisted approach to modeling 
the coarse-grained dynamics of such large, finite oscillator 
assemblies at and near full synchronization. The premise 
is that there exist a small number of coarse-grained vari- 
ables (observables) adequately describing the long-term 
dynamics, and that a closed evolution equation for these 
observables exists, but is not explicitly available. To ac- 
count for oscillator variability within the assembly, we 
treat both the variable oscillator properties (here, natu- 
ral frequencies ui) and the oscillator states (here, phase 
angles 0) as random variables. Recognizing a quick de- 
velopment of correlations between lo and 9, we express 
the latter as a polynomial expansion of the former (bor- 
rowing Wiener polynomial chaos (PC) tools [13); the PC 
expansion coefficients are our coarse observables. 

Availability of the governing equations for the variables 
of interest is a prerequisite to modeling and computation. 
We circumvent this step using the recently-developed 
equation-free (EF) framework for complex, multiscale 
systems modeling 153, In this framework we 

can perform system-level computational tasks without 
explicit knowledge of the coarse-grained equations; these 
unavailable equations are solved by designing, performing 
and processing the results of short bursts of appropriately 
initialized detailed (fine-scale, microscopic) simulations. 



We consider a paradigmatic model of coupled oscil- 
lators, the Kuramoto model, consisting of a population 
of A all-to-all, phase-coupled limit-cycle oscillators with 
i.i.d. to with distribution function g(u>). This model has 
been extensively studied because of its simplicity and cer- 
tain mathematical tractability, yet it is not merely a toy 
model. It appears as a normal form for general systems 
of coupled oscillators (e.g. Refs. jiot fTlj | ^ . 

We choose a Gaussian with standard deviation a u = 
0.1 for g(u>); however, our approach is not limited to this 
particular choice, nor to the Kuramoto model. Due to 
rotational symmetry, the mean frequency f2 = ^j^tOi/N 
can be set to without loss of generality. The governing 
equation for the phase angle of the ith oscillator 9i is 



d6 t K ^ 
— = ujj H > sir 



1< i < N, 



(1) 



3=1 



where K > is the coupling strength. Spontaneous syn- 
chronization (phase-locking) occurs at sufficiently large 
K. As K decreases across a critical value K tp , more and 
more oscillators desynchronize until they all essentially 
evolve with their own frequencies below another critical 
value K c 0, 0, 0] . Kuramoto introduced a complex 
order parameter re 1 ^ = SjLi e%6i to describe the long- 
time states; the effective radius r(t) measures the phase 
coherence; see also Ref. [22| for an order function. The 
asymptotic value of r (t — > oo) in the continuum-limit 
(A — > oo) exhibits a temporal analog of phase transition 
at K c 0. 

The order parameter r conveniently represents statis- 
tical behaviors around the critical point K = K c ; how- 
ever, r does not uniquely specify the microscopic state, 
and it may not adequately describe transient dynamics. 
The statistical moments of the phase angle distribution 

function M„ = (®3 ~ W zCili #i) > where n is 

a positive integer, are a "natural" first choice of coarse- 
grained observables (in a kinetic theory-like description) . 
Due to the symmetry, we consider only even-order mo- 
ments, and test whether a closure in terms of and 
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FIG. 1: Three coarsely identical (same M2 and M4) but 
microscopically different initializations (dashed, dotted and 
dot-dashed lines; see text) evolve along different trajectories, 
to a slow manifold and, ultimately, the same synchronized 
state (N — 300; K — 0.7 > Kt p ). Constraining the evolution 
to M2 = 0.017 (solid line) guides the trajectory directly to 
this slow manifold; the inset shows Mi becoming slaved to 
M2 during this constrained evolution by t « 2.0. 



FIG. 2: Time snapshots in the to — 6 plane for free evolu- 
tion (main panels; dashed line in Fig. and for constrained 
evolution (insets; solid line in Fig. respectively. Each dot 
represents an oscillator, and (a) to (d) are snapshots at t = 0, 
1.0, 2.0, and 6.0, respectively, marked by filled circles in Fig.Q 
Strong correlations develop during the initial transient stages 
( "oscillator sorting" ) . 



M.a is likely for K > K tp . We prepare several distinct 
initial phase angle distributions with identical coarse- 
grained values {M 2 = 0.017 and Ma = 0.0020); these 
phase angles are randomly assigned to oscillators. The 
phase portrait in Fig. ^ shows direct simulation [using 
Eq. for three of these initial assemblies; the trajec- 
tories are clearly distinct, suggesting that the dynamics 
cannot close simply on these two observables. Including 
higher order moments (such as M. 6 ) as observables does 
not remedy the situation. It is also clear, however, that 
the long-term dynamics lie on a low-dimensional man- 
ifold (ultimately a one-dimensional one) towards which 
all trajectories are eventually attracted. 

The dynamical differences among the three cases arise 
from the microscopic differences of the (macroscopically 
identical) initial conditions; this is best seen in the uj-9 
plane (Fig. [2J. Correlations between 9 and oj develop 
(the initial "cloud" in the lj-9 plane quickly evolves to 
a "curve"), as all transients initially approach the slow 
manifold: The oscillators "sort themselves out". These 
correlations were not accounted for when we assigned an- 
gles randomly to oscillators in the assembly. 

We now include a "remedial initialization" step, evolv- 
ing the dynamics by constraining them on prescribed val- 
ues of the moments, as a system of differential algebraic 
equations (DAEs) using Lagrange multipliers. The solid 
line in Fig. ^ shows this preparatory step with a con- 
straint on M.i only; constrained evolution brings the as- 
sembly down to the right point on the slow manifold, 
and the same "sorting" develops as in the aforementioned 
freely-evolving cases. Phase angle statistics alone do not, 
therefore, constitute good observables [2j|; to-9 correla- 
tions should be accounted for in the coarse description. 

Motivated by this observation, we explore the long- 



term dynamics with a different set of observables, treat- 
ing both 9 and lo as random variables. The former is 
expanded in Hermite polynomials of the latter, Gaus- 
sian random variable; Wiener polynomial chaos is the 
appropriate choice for Gaussian distributions. General- 
ized polynomial chaos (gPC) would be invoked for 
different frequency distributions (e.g. we also successfully 
used Legendre polynomial expansions for uniform g{uj)). 
For convenience, we introduce the normalized random 
variable £ = lu/ct u : 

i=0 i=0 * 11 *' 

where p is the highest order retained in the truncated 
series, (• , •) denotes the inner product with respect to 
the Gaussian measure, and Hi is the ith Hermite poly- 
nomial [Hq(x) = l,Hi(x) = x,H2(x) = x 2 — l,Hs(x) = 
x 3, — 3x, • • •]. Only odd-order on's are considered, due to 
symmetry. We will see that here the first two nonvan- 
ishing coefficients a% and 0:3 provide an adequate rep- 
resentation. Given a particular detailed realization of 
the oscillator state, its PC coefficients a^s are estimated 
through a least squares fitting algorithm, interpreting 9 
as an empirical function /(£) = ai£ + (£ 3 — 3£) and 

minimizing the residual R 2 = J2 3 ~ /(Cj '•'■> a ii a 3)] 2 ■ 
This procedure corresponds to the restriction (fine to 
coarse) step in the EF framework, described below. 

In the EF approach, appropriately initialized short 
bursts of detailed, fine-scale simulation are used to esti- 
mate quantities pertaining to the evolution of the coarse- 
grained variables (observables). Lacking an explicit 
coarse-grained model in terms of the first few PC coeffi- 
cients, we estimate the quantities necessary for scientific 
computation with it (time derivatives, action of Jaco- 
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FIG. 3: (color online) Coarse projective integration (dots) and detailed coupled oscillator dynamics (lines); N = 300. (a) Two 
PC coefficients (K = 0.4; full synchronization), (b) Two PC coefficients and a single "free" oscillator (K = 0.31). (c) Two 
PC coefficients and two "free" oscillators (K — 0.31). Natural frequencies are newly drawn from g(u>) at each lifting step (see 
text). Inset in (a): Schematic of a projective integration step: The last part (last two dots, at t n -i and t n ) of a short burst of 
direct integration (five dots) is used to estimate the local time derivative (solid line). PC values at a future time t = t n +M are 



"projected" through forward Euler, i.e., Qi| t= t 



tn). 



bians, residuals) through on demand numerical experi- 
mentation with the detailed, fine-scale model [Eq. JJ)]. 

The general procedure consists of (i) identifying good 
observables that sufficiently describe the coarse-grained 
dynamical state (here, a few c^'s), (ii) constructing a lift- 
ing operator, mapping the coarse description to one (or 
more, for variance reduction purposes) consistent fine- 
scale realization(s) [randomly drawing oj from g(w) and 
assigning 9, using Eq. © and given a; values], (Hi) evolv- 
ing the lifted, fine-scale initial conditions for certain time 
horizon, (iv) restricting the resulting fine-scale descrip- 
tion to the coarse observables [finding the PC coefficients 
of the final state], and (v) repeating the procedure as nec- 
essary to perform specific scientific computation steps. 
This is a general approach that has been combined with 
various fine-scale models 0,Hj|; see Refs. [iol l2l| . 

We first demonstrate coarse projective integration (2(|. 
Each group of five dots in Fig. represents the time hori- 
zon during which the detailed equations are integrated to 
enable the projective step; the local time derivatives of 
the observables are estimated here simply from the last 
two dots in each group. Coarse variable values at a pro- 
jected, future time are estimated using these derivatives 
and (for projective forward Euler) linear extrapolation 
in time [see the inset in Fig. 0(a)]. After the projection 
step we lift the coarse variables to consistent fine-scale 
realizations, and use these as the initial condition for an- 
other short burst of direct detailed integration [steps (ii) 
and (Hi) above]. Depending on the relative lengths of 
the projection step (MS) and the short run required to 
estimate the coarse time derivatives (nS), this procedure 
may significantly accelerate the computational evolution 
of the oscillator statistics; the cost of the lifting step (here 
negligible) must also be considered. At each lifting step, 
u> was newly drawn from g(u>), and the full integration 
(lines) and projective integration (dots) agree on the level 
of fluctuation among realizations. The PC coefficients 



display smooth behavior, nearly independently of partic- 
ular random draws; for the same random draw at every 
step, results would be even better. Projective integra- 
tion in Fig. 13 (a) reduced the computational effort in our 
illustrative direct integration by a factor of four. The nu- 
merical analysis of projective integration (stability, accu- 
racy, stepsize selection and estimation issues) is a topic 
of current research (see e.g., Refs. [13, E^, H3| ) ; here we 
simply demonstrated the procedure and its potential. 

Slightly below the transition value K tpi where only few 
oscillators become desynchronized, we consider the sys- 
tem as a combination of synchronized "bulk" and a few 
"free" oscillators. Good coarse-grained observables then 
are a few PC coefficients for the "bulk" synchronized os- 
cillators and the phase angle(s) of the (few) desynchro- 
nized one(s). The EF approach can be directly "wrapped 
around" this alternative representation. Both for the 
one free and two free oscillator cases, projective inte- 
grations on the new observables successfully track (and 
accelerate) direct detailed simulations [Figs. (b) and 
(c)]. These "good observables" are suggested by direct in- 
spection and common sense; for more complicated, high- 
dimensional problems, good state parameterizations re- 
quire modern data mining algorithms. Diffusion maps on 
graphs constructed by the data [2^] are a promising tool 
for detecting good "reduction coordinates" (observables) 
on which to base EF computations. 

Direct, long-time simulation is often inefficient in com- 
puting long-time (stationary) states. Numerical bifurca- 
tion algorithms, more appropriate for stability and para- 
metric analysis, can be implemented in an equation-free 
framework: The residual and the action of the unavail- 
able Jacobian are numerically estimated through short 
bursts of appropriately initialized detailed simulations. 
Starting from a coarse-grained initial condition, we lift, 
and integrate the full model for time AT. We then re- 
strict to the observables of the final state $at; this is the 
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FIG. 4: (color online) Coarse bifurcation diagram for the full 
synchronization regime (K > Ktp), obtained through coarse 
Newton-GMRES method and pseudo arc-length continuation 
(N = 300). The same variables as in Fig. [3] (b) are used; 
the phase angle of the single "free" oscillator (9^ ree ) is an ex- 
tra observable (its natural frequency is positive in this case). 
The PC coefficients, obtained by discounting the "free" os- 
cillator, exhibit nearly the same values both for the stable 
(filled circles) and the unstable (open circles) branch (only a± 
is shown here). Only 0f ree shows significant variation along 
the two branches. Arrows are included to guide the eye. 



coarse time-stepper. We now solve for the fixed point sat- 



isfying $at 




0, using the coarse 



Qjree j \ Qfree 

Newton-GMRES |28(, a matrix-free iterative method (to- 
gether with the pseudo arc- length continuation); addi- 
tional coarse observables 9^ ree are appended when nec- 
essary. We construct bifurcation diagrams like the one in 
(Fig. 01 with respect to the parameter K, showing a turn- 
ing point (actually, a "sniper") bifurcation at K = K tp . 
A single oscillator (whose phase angle Qf ree is treated as 
a separate coarse observable) becoming "free" from the 
synchronized "bulk" at that point. For sufficiently large 
K values (when r « 1) analytical estimates of certain 
elements of the shape of the u>-6 correlation become pos- 
sible (e.g., from Eq. QJ one can obtain, at steady state, 
ai w a^/K, in reasonable agreement with our steady 
state computations at K >~ 0.5). 

In summary, the EF multiscale approach was success- 
fully used for coarse-grained dynamic computations of 
finite assemblies of non-identical coupled oscillators; the 
derivation of explicit closures at- and close to the syn- 
chronization regime was circumvented. Initial transient 
"sorting" of the oscillators, establishing correlations be- 
tween natural frequencies and phase angles, suggested 
Wiener PC coefficients as the appropriate coarse observ- 
ables. If the problem dynamics can be coarse-grained, 
traditional numerical analysis algorithms can be used as 
protocols for the "intelligent" design of short bursts of 
computational experiments with the detailed, fine-scale 
model. The approach can be directly generalized to an- 
alyze the simulation and modeling of more complicated 
oscillator dynamics. 
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